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I. Introduction 


This project involves obtaining GPS measurements in Scandinavia, and using the measurements 
to correct tide-gauge measurements for the rebound effect, and to estimate the viscosity profile of the 
Earth's mantle. Below, we report on several aspects of this project. 

II. GPS Measurements 

The BIFROST permanent networks set up by Onsala Space Observatory and the Finnish Geodetic 
Institute continues to operate, and the data are continuously being analyzed. 

The planned analysis of the BIFROST GPS data was carried out. In March, we produced a new 
velocity solution. 

III. Geodetic and Geophysical Results 

These have been reported in a recent submission to JGR, Continuous GPS measurements of post- 
glacial adjustment in Fennoscandia, 1. Geodetic results. This manuscript is attached as Appendix A. 

IV. Sea Level and Postglacial Rebound WWW Site 

This site has been down since May 30 for extendive remodeling. We are implementing our post- 
glacial rebound calculator with several hundred possible model combinations. 



Appendix A 


Attached is the paper submitted to the Journal of Geophysical Research 
GPS measurements of postglacial adjustment in Fennoscandia, 1. Geodetic 


in April, 2000: Continuous 
results. 
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Abstract. Project BIFROST (Baseline Inferences for Fennoscandian Rebound 
Observations, Sea-level, and Tectonics) combines networks of continuously operating 
GPS receivers in Sweden and Finland to measure ongoing crustal deformation due to 
glacial isostatic adjustment (GIA). We present an analysis of data collected in the 
years 1993-1998. We compare the GPS determinations of three-dimensional crustal 
motion to predictions calculated using the high resolution Fennoscandian deglaciation 
model recently proposed by Lambeck et al. [ 1998a, b]. We find that the the maximum 
observed uplift rate (~10 mm yr -1 ) and the maximum predicted uplift rate agree to 
better than 1 mm yr -1 . The patterns of uplift also agree quite well, although differences 
are discernible. The \ 2 difference between predicted and GPS-observed radial rates is 
reduced by a factor of 5-6 compared to that for the “null” (no uplift) model, depending 
on whether a mean difference is first removed. The north components of velocity agree 
at. about the same relative level, whereas the agreement for the east, components is 
worse, a problem possibly related to the lack of bias fixing. We have also compared 
the values for the observed radial deformation rates to those based on sea-level rates 
from Baltic tide gauges. The weighted RMS difference between GPS and tide-gauge 
rates (after removing a mean) is 0.6 mm yr -1 , giving an indication of the combined 
accuracy of the GPS and tide-gauge measurement systems. Spectral analysis of the 
time series of position estimates yields spectral indices in the range —1 to —2. An EOF 
analysis indicates, however, that much of this power is correlated among the sites. The 
correlation appears to be regional and falls off only slightly with distance. Some of this 
correlated noise is associated with snow accumulation on the antennas or, for those 
antennas with radomes, on the radomes. This problem has caused us to modify the 
radomes used several times, leading to one of our more significant sources of uncertainty. 
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1. Introduction 

The last 800 ka of the current ice age have been characterized by a series of ''glacial 
cycles.” each with a period of approximately 100 kyr [e.g.. Broeker and van Donk. 
1971]. Within each cycle a relatively slow glaciation phase, culminating in massive 
ice complexes over most of the high latitude continental regions, was followed by a 
much more rapid deglaciation event. For example, during the last glacial maximum, 
which occurred just ~20 kyr B.P., the ice sheets reached thicknesses of 2-3 km or 
more in Fennoscandia, Canada, Antarctica, Greenland, Siberia, and Arctic Canada 
[Denton and Hughes. 1981]. Remarkably, in a matter of only 15 kyr, a large proportion 
of ice disintegrated, leaving most of these regions ice-free and raising world-wide 
ocean levels by over 100 m [e.g., Chappell and Shackleton , 1986]. The redistribution of 
surface ice-water mass implied by these glaciation/deglaciation episodes has induced an 
appreciable and ongoing isostatic adjustment of the planet. 

Glacial isostatic adjustment (G1A) is manifested in a wide variety of past- and 
present-day geophysical observables that have previously been used to study this 
process, including time series of ancient sea-level elevations (relative to present-day sea 
level), the modern tide gauge record, gravity anomalies, present-day secular variations in 
the global gravity field, and present-day secular variations in the Earth's rotational state 
(variations in length of day and motion of the rotation pole) [e.g., Milne , 1998]. These 
observations provide an indirect inference of present-day ongoing crustal deformation. 
Direct high-accuracy measurements of crustal deformation, even in a particular region, 
were not possible, however, before the advent of space geodetic techniques, and even 
these techniques have only quite recently been able to achieve the required accuracy. 
Very long baseline interferometry (VLBI) has been available as a high-accuracy geodetic 
technique for over 20 years; several studies have now addressed the effects of GIA on 
VLBI determinations of site velocities [e.g., James and Lambert , 1993; Mitrovica et al. , 
1993]. Unfortunately, the global VLBI network is extremely sparse, and only 2-3 sites 
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exhibit much sensitivity to GIA. A small error in the velocity determination of these 
sites will unduly influence conclusions regarding mantle viscosity and ice-sheet histories 
[Mitrovica et al., 1993; Mitrovica et al., 1994b]. 

Geodesy with the Global Positioning System (GPS) affords us several advantages 
relative to geodesy with VLBI related to the relatively low cost of the GPS receivers 
(~$15 K or less). As a result of this low cost it is feasible to deploy a dense network 
of receivers across a region. The detailed pattern of deformation may thus be inferred. 
The low cost of GPS receivers moreover makes it financially feasible to dedicate a 
group of GPS receivers to ‘'permanent” sites within a region in order to acquire 
“continuous” measurements. In principle, given sources of error that axe sufficiently 
steady state, it should be possible to “beat down” the noise and thus determine estimates 
of velocity much more quickly and accurately than with conventional “campaign” 

GPS measurements. The number of GPS receivers in permanent GPS networks is 
quickly outpacing the number used for conventional campaign high-accuracy geodetic 
measurements [e.g., Segall and Davis, 1997]. 

Geodesy with GPS has been steadily and significantly improving in precision and 
accuracy over the past ten years. This improvement has mainly been due to advances 
in the GPS satellite constellation, GPS receiver design, and analysis techniques. The 
demonstrated repeatability of horizontal position estimates obtained from GPS is 
currently at the few-mm level on regional and local scales and at the 10-mm level on 
global scales [e.g., Blewitt , 1993]. The repeatability in the vertical baseline component 
is typically 3-5 times worse. The level of accuracy achievable in a single day with the 
GPS technique is thus in principle equal to that achievable with VLBI. 

In August 1993, we established a network of permanently operating GPS receivers 
in Sweden [BIFROST Project , 1996]. A number of sites in Finland were also temporarily 
occupied. In 1994-6, permanent GPS sites were established in Finland. The GPS sites 
within these networks, along with several already existing sites of the International GPS 
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Network for Geodynamics (IGS) sites in Norway operated by the Norwegian Mapping 
Authority, make up a dense regional Fennoscandian GPS network of an intersite spacing 
of ~100 km and a total area of ~2000 x 2000 km, situated within the area covered 
by ice at the last glacial maximum. Investigators from three Nordic and two North 
American institutions formed Project BIFROST (Baseline Inferences for Fennoscandian 
Rebound, Sea level, and Tectonics) [BIFROST Project, 1996]. One of the primary goals 
of BIFROST is to use the three-dimensional velocity vectors from the BIFROST GPS 
network to provide a new GIA observable for the determination of Earth structure and 
Fennoscandian ice history. In this paper, we report on the first results from this effort. 
\Ve will present the BIFROST GPS networks and data sets, and describe the analysis 
of the GPS data. We include a thorough discussion of errors since the GIA signals we 
are attempting to measure are quite small (sub-cm). 

2. The BIFROST GPS Networks 

The BIFROST GPS networks (Figure 1) are composed of the permanent GPS 
networks of Sweden (SWEPOS™) and Finland (FinnRef). Table 1 describes the 
histories of the BIFROST sites. This table contains the dates and configurations of the 
original installations, as well as the dates on which significant, modifications were made 
to the site hardware. (Minor changes, such as those affecting communication only, are 
not indicated in Table 1.) The hardware described in Table 1 includes the GPS receiver 
type, GPS antenna type, radome type, monument type, and approximate positions. 
IERS Domes Numbers for those sites which are IGS sites are given as well. 

Below, we discuss the individual aspects of both BIFROST networks. 

2.1. The SWEPOS GPS Network 

The Swedish nationwide multipurpose network of twenty-four permanent GPS 
stations, SWEPOS, was established in 1993. On July 1, 1998. SWEPOS attained 
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full operational capability for real-time positioning at the meter accuracy and for 
post- processing applications with centimeter accuracy. Real-time positioning at the 
centimeter/decimeter level is planned for 2002. The National Land Survey of Sweden 
(LMV) is responsible for the maintenance and the operation of the SWEPOS network 

The SWEPOS network (Figure 1) currently consists of 21 continuously operating 
GPS stations. The “standard" SWEPOS monument (designed by the LMV and denoted 
as “SWEPOS” in Table 1) consists of a 3-m tall concrete circular pillar atop a concrete 
platform. At five sites (Kiruna, Lovo, Martsbo, Norkoping, and Skelleftea) a second 
pillar is available to serve as an alternate and as a platform for test measurements. The 
pillars are built on bedrock and the line of sight from the top to the GPS satellites is 
unblocked at elevation angles above 10° and often lower. Each pillar is supported by four 
internal steel rods set 1 m into the underlying rock. Heating coils are helically around 
each concrete pillar. Insulating material consisting of helically wound corrugated plastic 
sheet and rockwool surrounds the wire-wrapped pillar. A temperature sensor is fit into 
a small cavity inside the pillar and is connected to a thermostat unit in the instrument, 
cabin. The thermostat maintains the temperature of the sensor above 15 °C. On the top 
of each pillar is a plate for the attachment of the GPS antenna, tribrach, and adaptor. 

The pillars at each site are surrounded by a network of steel pins, driven into the 
rock so that their tops protrude a few centimeters above the surface. This local network, 
covering an area of approximately 15 m x 15 m, is used to monitor the stability of 
the concrete pillars. The GPS antenna is removed from the pillar and replaced with a 
theodolite, which is used to measure the horizontal and vertical angles to the steel pins. 
Through resection the position of the pillar can be calculated. In this manner, the local 
position and orientation of the pillar may be monitored to better than 1 mm. The first 
such measurements were obtained during summer 1993 and repeated annually except 
for the monument in Leksand where measurements are carried out monthly. The results 
of these measurements are described below. 



Several sites have slight variations in monumentation. The Onsala site has a 
different monument, due to its earlier construction as an IGS site. The Onsala 
monument consists of a 1 m tall pillar with a square cross-section and without heating 
control or insulating material. The Jonkoping pillar is 1 m shorter than the standard 
pillar (for air traffic safety), and also is not heated. The Lovo and Martsbo monuments 
are SWEPOS monuments built over pre-existing pillars of rectangular cross-section. In 
addition, the multipath environment at these two station might be worse than at others 
due to pre-existing construction. 

Each SWEPOS site has a 3 m x 2 m hut housing the GPS receivers, backup 
batteries, computer, and Internet connection. All SWEPOS sites are equipped with two 
or more GPS receiver systems, and AOA SNR-8000 and an Ashtec Z-XII, connected to 
a single Dome- Margolin type antenna. At four stations dual-frequency GPS/GLONASS 
receivers are also installed. The stations are equipped with a power backup system, 
which can run the station for 48 hours if the main power fails. All stations are connected 
to the control center via leased 64 kbit lines and a redundant 19.2 kbit X.25 line. 

2.2. The FinnRef GPS Network 

Planning for the FinnRef network (Figure 1) started at the Finnish Geodetic 
Institute (FGI) at the end of 1992, when it was decided that a network of 12 stations 
would be established. Possible site candidates were chosen with several criteria in mind: 
(1) the network should cover the country so that the maximum land uplift differences 
could be sampled; (2) the stations should be built on bedrock and there should be open 
sky above an elevation angle of 15°; (3) absolute gravity has been or can be measured 
on the spot; and (4) stations should easily be connected to the precise levelling network 
and to the telephone and electricity networks. The criterion (2) has generally, but 
not universally, been met, and in most cases the horizon is 10° or lower. Planning, 
construction and use of the network are described in more detail elsewhere [e.g., Koivula 
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etal . , 1997, 1998] 

At the Joensuu, Kuusamo, Vaasa, Virolahti, Olkiluoto, Kivetty and Romuvaara 
stations we constructed heated wooden cabins of area 1.5-2 m x 3 m to house the GPS 
and other electronics. Existing buildings were used at all other sites. 

Some stations are located close to other institutions where on-site personnel can 
assist in case of minor problems. Kevo is on the premises of the Subarctic Research 
Center of the University of Turku, Oulu is at the Aarne Karjalainen Observatory of 
the University of Oulu; Sodankyla is visited weekly by local staff of the Sodankyla 
Geophysical Observatory; and Tuorla is at the Astronomical Observatory of the 
University of Turku. At Kivetty, Olkiluoto and Romuvaara there are also contact 
persons who can check the stations. Metsahovi is at the Space Geodetic Observatory of 
the FGI. 

Three different types of antenna platforms are used for FinnRef. The standard 
configuration is a 2.5 m high steel grid mast, which is used at Joensuu, Kuusamo. 
Sodankyla, Tuorla, Vaasa and Virolahti. A similar mast, but 5 m high is used at Kevo. 
In the case of the 2.5 m mast, the thermal expansion effects amount to a height variation 
of less than ±1 mm during the annual temperature cycle. This variation is considered 
to be acceptable. Around the antenna masts, there are three reference bench marks, 
and connections to the first order levelling network have also been established. 

Two stations have higher masts. There is an anchored 25 m high steel grid mast at 
Metsahovi and at Oulu there is a cylindrical 8 m steel mast. In both cases the height of 
the GPS antenna is stabilized with an invar rod. The antenna is isolated from the mast 
with an attachment piece and a spring system, which is anchored to the bedrock with 
an invar rod or wire [Paunonen, 1993]. The system for Oulu was adapted from that at 
Metsahovi. Three stations, Olkiluoto, Kivetty and Romuvaara, were built in cooperation 
with Posiva Oy, a company which is responsible for locating sites for disposal of nuclear 
waste. Local networks around these sites are remeasured semiannually in order to locate 
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possible deformations. For this reason more stable concrete pillars were chosen at these 
sites [ Chen and Kakkuri , 1994]. 

All stations are equipped with Ashtech Z-XII GPS receivers, Dorne-Margolin type 
antennas, modems and power supplies. The exception is Metsahovi where an AOA 
SNR-8100 receiver is in use. At Metsahovi there is also an external H-maser; at all 
other stations the receiver’s internal oscillator is used. Except Metsahovi and Tuorla. 
all antennas are equipped with a radome. These have proven less than satisfactory for 
their stated purpose, but it was decided not to change the antenna mount further due 
to the experience with SWEPOS. 

Data are collected using a sampling interval of 30 s and a 5° elevation-angle 
cut-off. During the 1998/9 winter CB00 software was installed into all receivers and 
the sites were equipped with Vaisala PTU 220 meteosensors. The station histories are 
summarized in Table 1. 

3. Data Analysis and Geodetic Results 

The dual-frequency GPS phase and pseudorange data were processed using the 
2nd release of GIPSY software developed at Jet Propulsion Laboratory (JPL) [e.g., 
Webb and Zumberge, 1993]. Dual-frequency phase and pseudorange data from a single 
30-hour period acquired from all the sites in the network are analyzed simultaneously. 
(There is a 3-hour overlap at each end of each observing session.) The GPS data 
are decimated to achieve an effective sample rate of 300 sec; decimation is performed 
to maintain a manageable level of utilized disk space. For each 30-hour data set we 
estimated the usual set of parameters, including oscillator (‘Flock”) corrections, site 
positions, atmospheric zenith delay parameters, and ambiguity parameters. Satellite 
orbit parameters were highly constrained to the values distributed by the IGS based on 
a solution involving a global network of GPS sites. Temporal variations in the clock and 
atmosphere parameters are modeled as independent random walks [Webb and Zumberge , 
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1993]. We adopted a minimum elevation angle of 15° for all stations. Corrections for 
the motion associated with ocean loading and solid-Earth tides were incorporated in the 
model. 

For the SWEPOS sites, which now have multiple antennas, the SNR-8000 data 
are used in the solution up until August 1, 1998, and the Z-XII data thereafter. For a 
period following this date we performed a number of solutions with both GPS receivers 
and determined differences at the 1 mm level or less. 

We adopted the value of 10 mm for the uncertainties in the phase measurements at 
each frequency. The instrumental uncertainties for such measurements are much smaller, 
perhaps 1-3 mm [ Spilker , 1996]. However, experience within the GPS community has 
shown that the scatter of the time series is greater than the theoretical value based on 
instrumental noise only. The increase in the scatter above the predicted value can of 
course be attributed to unmodeled phase variations, which may or may not have a white 
noise (or nearly white noise) nature. In Section 5 we discuss a number of errors which 
might contribute to this increased scatter, and we investigate the spectral characteristics 
of the site position variations. It is important to remember throughout the paper, 
though, that the uncertainties for the estimated parameters, including site position and 
therefore velocity and geophysical parameters, are approximations. 

Data processing utilizes a ‘‘no-fiducial” technique described by Heflin et al. [1992] 
wherein station coordinates have weak a priori constraints. The results presented in the 
following sections were achieved without fixing the estimated phase biases to integer 
values, since the software cannot handle automatically such an extensive data set. 
Independent tests with a smaller data set (SWEPOS stations only) have shown that 
“bias fixing” leads to a decrease in the formal errors of about 20-30%. 

All geodetic positions obtained in the GIPSY analysis are finally transformed into 
the International Terrestrial Reference Frame of 1996 (ITRF96) [Sillard et al. , 1998], 
This transformation creates a slight inconsistency since the satellite orbits prior to 
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September 1996 are referred to different editions of ITRF (ITRF92 until July 1994, 
ITRF93N until June 1995, and ITRF94 until September 1996). The problem that 
concerns us here is that the frames have slightly different net rotations and translations. 
Using a set of core stations constituted by those that are jointly present in pairwise 
successive reference frames and applying weights as given by the velocity uncertainties, 
the rotations and translations derived by least-squares adjustment can amount to 
1 mm yr~ 1 at the Earth’s surface in any component. We have therefore estimated and 
applied the inter-frame rotation and translation parameters to correct the time series of 
station positions for these biases. The results are thus determined in a rigid frame that 
is co-moving with the ITRF96 frame, which for our sites is dominated by the ITRF96 
realization of the Eurasian plate motion. We, however, are interested in deformations 
relative to the Eurasian plate motion. 

Removing the motion of the Eurasian plate requires a decision regarding the type 
of motion the co-moving Eurasian frame should be allowed. The requirement to observe 
deformation from a rigid co-travelling frame implies the need to suppress a deformation 
mode conveyed by a scale rate parameter. Although any residual rotation in one corner 
of the region could easily be corrected for by adding small rigid rotations, the case is 
slightly more intricate when one simultaneously considers translations. Considering 
horizontal motion, the virtue of the GPS data in application to GIA lies in the ability to 
resolve intersite motions; subtracting the motion of a rigid frame will have no influence 
on relative deformation. In the case of radial motions, GPS data will not only be used to 
study relative deformation, but also offers the prospect of studying the vertical motion 
of the crust, for example, in comparison with the sea surface. An "absolute” frame is 
therefore desired. 

Allowing for translations in the co-moving frame will effect the estimates of vertical 
components of site motion. One could argue that the European stations are well 
established and stable, so that the ITRF96 motion after correction for vertical motion 
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from one or a set of models would provide a proper regional vertical reference. The 
associated co-moving frame would be constructed by estimating rotation and translation 
rate parameters using the subset of velocity estimates at those European IGS stations 
that also were used in the projection stage. 

On the other hand, the large number of stations world-wide now included in the 
ITRF should provide a stable constraint. In the larger, global set, local vertical motion 
will appear as less correlated. Accepting this argument, the consequence is to not allow 
relative vertical motion between the co-moving frame and the ITRF, i.e., use only the 
horizontal projection of the ITRF site motion vectors. This method brings about an 
advantage, namely that the observed data do not have to be “corrected” with a model, 
thereby avoiding problems of circular arguments at the stage where the network rates 
are interpreted. The associated co-moving frame is simply constructed by solving for 
rotation rates only. We have therefore adopted this method. 

The analysis described above yields a time series for each station of three- 
dimensional position in the ITRF96 “Eurasia co-moving” reference frame. Time series 
for the BIFROST stations and Tromso, a nearby IGS site, are shown in Figure 2. 

3.1. Antenna-Related Issues 

Several of the time series in Figure 2 for sites of the SWEPOS network exhibit 
one or more “jumps.” We do not believe that these jumps represent motions of the 
GPS antenna. The jumps are associated with (1) removal and replacement of the GPS 
antennas, (2) changes in antenna radomes, and (3) rapid changes in snow accumulation. 
Snow accumulation is discussed in Section 5.4. 

In the removal and replacement of the GPS antennas to perform the local site 
surveys mentioned above, the GPS antenna is positioned on the monument by means 
of a threaded hole and a standard 5/8” surveyor’s bolt attached to a metal plate that 
has been permanently set into the concrete at the top of the pillar. The GPS antenna 
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is screwed onto the bolt until it refuses to rotate. When the antennas are removed 
and replaced, the orientation of the antenna is checked to insure that it is the same as 
before removal. Given that the surveyor’s bolt has 5 threads per inch, a rather large 
orientation error of 45° would lead to a vertical displacement of only 0.6 mm, and to no 
horizontal displacement. The “jumps” in Figure 2 associated with these surveys, on the 
other hand, can be at the 10 mm level. 

A likely explanation for these jumps is that very small differences in antenna 
orientation lead to changes in phase errors because of electromagnetic coupling [Eldsegui 
et al . , 1995; Jaldehag et al . , 1996b] and antenna phase center variations [ Schupler et 
al . , 1994]. Both these sources of error are potentially elevation- and azimuth-angle 
dependent, and in the case of the former the position relative to the pillar and metal 
plate is critical. If one imagines that the phase errors induced by these two phenomena 
can be represented as a series of spherical harmonics, with angular arguments of azimuth 
and elevation angles, then the contribution from the £ — 1 term is indistinguishable 
from the contribution to phase variations from a site position offset. 

Other apparent jumps occur when there were changes in the antenna radomes. The 
original radomes installed on the SWEPOS sites were designed at Delft University of 
Technology. During the winter of the first year of the experiment, snow accumulated 
significantly on these radomes, and our observations led us to believe that this 
accumulation could be reduced by a redesigned radome having no horizontal surfaces. 
Re-designed radomes (“Type A”) were emplaced in the winter and spring of 1995. We 
later discovered that the paint process used on these radomes were defective. These 
radomes were thus removed in the spring and summer of 1996 and later that year 
replaced with improved radomes (“Type B”). Each of these changes appears to produce 
offsets in the time series. The radomes are discussed in detail by Emardson et al. [2000] . 

As an ad hoc treatment for these errors, we have simply estimated three-dimensional 
offsets in position on the epochs at which radomes were changed, the GPS antennas 
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removed and replaced, or the antenna rotated. These changes are summarized in 
Table 1. The site velocity was assumed to be constant for the entire experiment. This 
ad hoc procedure is not very satisfying, since the existence of the offsets is an indication 
of an error source which could conceivably have a temporal variation and therefore could 
effect the estimate of the rate. 

3.2. Determination of Station Velocities 

In this section, we report and compare several methods for determining the station 
velocities. Since, after the analysis described above, the time series of station positions 
are in a consistent reference frame, it is in principle simply a matter of fitting a straight 
line component by component and site by site to the time series. This method does not 
yield determinations of the correlations of the errors in the rate estimates, but these are 
formally very small since the orbit parameters were highly constrained in the original 
solutions. In order to gain a quantitative understanding of the effects of errors that 
are difficult or impossible to model, we present several different analyses for the rates. 
Each analysis uses different models for the variation of site position with time as well as 
different editing criteria. 

In the following, the standard deviations we report are the so-called “standard 
errors.” These standard errors are based on the phase uncertainties used in the daily 
least-squares analysis, described above, propagated through that analysis to yield 
standard deviations for daily determinations. (The standard errors are the uncertainties 
shown in Figure 2.) The errors for estimates obtained on different days are assumed to 
be uncorrelated. These standard deviations are then propagated through the analyses 
described below to yield the standard errors for the rate parameters. In general, the 
reduced x 2 postfit residuals are close to unity, indicating a reasonable fit, but this 
statistic may not be an accurate measure of the accuracy of the rate estimates. In the 
next section, we assess these standard errors and the accuracy of the rate estimates. 
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The solutions for the three-dimensional crustal deformation velocities and their 
standard errors for the different analyses are presented in Table 2. In all cases data 
prior to 1998.0 only were used. Below, we describe the different solutions. 

3.2.1. Standard solution. In the “standard solution,'’ sites having a short time 
series are analyzed differently from those having a longer time series. For longer time 
series (generally those with a timespan of two years or greater), the model for the 
position estimates included: a mean value, a constant rate, an admittance parameter 
for atmospheric loading [vanDam and Wahr , 1993], and periodic terms with frequencies 
of one, two, and three cycles per year. (The periodic terms are meant to model 
approximately the effects of snow accumulation.) For the short time series, no periodic 
terms were included. For the radial (“up”) components- only, offsets at each antenna 
change (Table 1) were included. No editing was performed. 

3.2.2. Edited solution. The parameterization is identical to the standard analysis, 
with the following exceptions. No difference was drawn between short and long time 
series. No admittance parameters for atmospheric loading were estimated, and the 
annual periodic, term only was included. Offsets for antenna changes were estimated 
for all components. An editing loop was included that deleted data whose postfit 
residual was greater than three times the weighted root-mean-square (WRMS) residual. 
This loop was repeated three times. Data from time periods 1995.000-1995.104, 
1995.370-1995.520, and 1996.438-1996.616 were automatically deleted. These periods 
represent timespans during which radome change/replacement was occurring across the 
network. 

3.2.3. EOF solution. An empirical orthogonal function (EOF) analysis for the 
radial rates was performed on a subset of the data. This analysis is described fully in 
Section 5.1. Table 2 includes the radial rates resulting from this solution. 
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3.3. Comparison of Solutions 

In Figure 3 we have compared rate determinations for all three components for the 
Standard and Edited solutions, and the radial rates for the EOF solution. In general the 
agreement is quite good despite the differences between these analyses, but in several 
cases the disagreement is quite large. Generally, the extreme differences occur for rates 
with larger error bars, indicating a shorter time series. Thus, these rates are relatively 
less stable to changes in analysis strategy. An obvious and important exception to 
this observation is the rate determination for Kiruna. From Table 2, we see that the 
radial rate for this site from the Standard solution is 11.7 ± 0.9 mm yr -1 , whereas 
the rate from the EOF solution is 10.6 ± 1.9 mm yr -1 and the rate from from Edited 
solution is 4.7 ± 0.9 mm yr -1 . Kiruna, a site in the north of Sweden, has a great deal 
of snowfall. The average winter precipitation is over 100 mm with an average winter 
temperature of approximately —10 °C [Schonwiese and Rapp, 1997]. Kiruna was the site 
at which snow accumulation was first noticed, and it has some of the largest seasonal 
signatures (Figure 4). Furthermore, only in the latest period has Kiruna had more than 
one continuous year of data with no offsets (Figure 2). Thus, Kiruna may be the most 
outstanding example of the types of problems discussed above. 

The differences between the Standard and Edited solutions have WRMS differences 
of 0.9 mm yr -1 (east), 0.5 mm yr -1 (north), and 1.6 mm yr -1 (radial). Treated as 
independent solutions of equal weight, this would indicate a typical uncertainty for a 
single rate of 0.6 mm yr -1 (east), 0.4 mm vr _1 (north), and 1.1 mm yr -1 (radial). This 
represents a scaling of the typical standard errors of 3-4 for the horizontal components 
and 1.5 for the radial. However, these solutions are not independent and in fact use 
much the same data set. From this point of view, the differences between the standard 
and edited solutions are large. 

The accuracy of the rates estimates cannot be perfectly assessed by a comparison 
like that of Figure 3, since the analyses share overlapping data sets and thus the 
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comparisons will not reveal common errors. Nevertheless, Figure 3 gives us an indication 
of the shortcomings in our modeling of the time series. By comparing to models 
for crustal deformation, we can obtain a rough assessment of our errors, and this 
comparison is contained in Section 4. Sources of error other than those discussed above 
are analyzed in Section 5. Phenomena that may actually influence the position of the 
site are considered also in Section 4. A more complete understanding of our errors will 
come with longer time series having no equipment changes are established. 

4. Interpretation of Observed Deformation Rates 

As stated above, a primary goal of the BIFROST Project is to provide a new 
and useful GIA observable with which to constrain models of the GIA process in 
Fennoscandia. In order to achieve this goal, the observations must exhibit a coherent 
signal that is clearly related to the regional GIA process. In this section we test this 
requirement by comparing the observed three-dimensional deformation rate signal to 
numerical predictions of this field and to the apparent sea-level signal that has long 
been associated with the GIA process. We also consider several other geophysical effects 
that may produce temporal variations in site position. 

4.1. Glacial Isostatic Adjustment 

A number of publications, some dating back to the 1930 : s [e.g., Haskell , 1935; 
Vening Meinesz, 1937], have employed sea-level observations to infer Earth viscosity 
and ice sheet parameters in the Fennoscandian region [e.g., Fjeldskaar, 1994; Mitrovica . 
1996; Lambeck et al. 1998a; Davis et al 1999]. In the recent study of Lambeck et al. 
[1998a], a three-layer Earth viscosity model and a regional ice model were proposed that 
provide a good fit to a carefully compiled and extensive data set based on geological 
sea-level markers. The preferred Earth models are defined by a. lithospheric thickness of 
65-85 km, an upper mantle viscosity of 3-4 x 10 20 Pa s and a lower mantle viscosity 
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that is a factor of ten or more greater than the upper mantle value. This range of 
three-layer Earth models and ice model were also found to produce a good fit to recent 
instrumented sea- level and lake-level records [Lambeck et al., 1998b]. (This more recent, 
shorter time-scale data apparently did not allow a robust inference of lower mantle 
viscosity.) 

A number of previous inferences that appear to disagree with the above described 
viscosity profile [e.g., Wolf. 1987; Fjeldskaar, 1994] are, in fact, found to be compatible 
when the resolving depth of the various data sets is considered [Mitrovica, 1996]. The 
inference of Mitrovica and Peltier [1993], which is based on the so-called Fennoscandian 
relaxation spectrum [McConnell, 1968], is not consistent with the Lambeck et al. [1998a] 
result. However, recent studies show that the paleoshoreline data upon which this 
spectrum is based require some revision [Wolf, 1996]. Indeed, a recent re-analysis of the 
spectrum was found to eliminate the inconsistency [Wieczerkowski et al., 1999] between 
these two inferences. We have limited the above discussion to recent GIA analyses that 
considered data from north western Europe. We have not considered recent analyses 
based on global sea-level data sets (which may include data from north western Europe) 
in order to avoid the potential bias introduced to these inferences from lateral variations 
in viscosity structure. 

The following predictions are based on a spherically symmetric, compressible, 
Maxwell viscoelastic Earth model. We choose a three-layer viscosity model defined 
by a lithospheric thickness of 70 km, an upper mantle of 4 x 10 20 Pa s and a lower 
mantle viscosity of 5 x 10 21 Pa s. (These values are consistent with the sea-level 
constraints discussed above.) The elastic structure of our Earth model is taken directly 
from the seismically constrained PREM [Dziewonski and Anderson, 1981]. The ocean 
component of the surface load is computed via a revised sea-level algorithm that 
solves the sea-level equation [Farrell and Clark, 1976] in a gravitationally self consistent 
manner while incorporating the effects of GIA-induced perturbations to the Earth’s 
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rotation vector [e.g., Milne , 1998] and the postglacial influx of ocean water/ meltwater to 
once-ice-covered regions [Milne, 1998]. Predictions of the load-induced three-dimensional 
deformation rate signal are calculated via the theory of Mitrovica et al. [1994a]. This 
theory has also been extended to incorporate the influence of GIA-induced perturbations 
in the Earth’s rotation vector Mitrovica et al. [2000]. The relative importance of the 
different components of the model will be described in a future publication. 

We require an ice model that provides a good fit to the sea-level observations for 
our choice of Earth model. This criterion is met by the model proposed by Lambeck et 
al. [1998a]. However, in order to accurately solve the sea-level equation and realistically 
compute GIA-induced perturbations to the Earth’s rotation vector, we require a global 
ice model. To meet both of these requirements we remove the Fennoscandian and 
Barents Sea components of the lower resolution, global ICE-3G [ Tushingham and Peltier. 
1991] model and replace these by the high resolution, regional model proposed by 
Lambeck et al. [1998a]. The contours of uplift that we calculate using this ice model and 
the Earth model described above are shown in Figure 4a. 

In order to illustrate the pattern of uplift that we observe from the BIFROST 
network, we have fit a simple surface to the vertical rates from the standard solution. 
The model we have chosen for the radial rate it at longitude A and latitude 0 is a 
two-dimensional Gaussian model: 

it( A, 0) = m 0 + ii a exp [— (A — A 0 ) 2 -I- u' 2 (A — A o )(0 — 0 O ) + (0 — 0 O ) 2 }] (1) 

The Gaussian is centered at (A 0 , 0 O ) and has maximal value ii 0 -I- it a and minimal value 
ii 0 . The parameters w\ and u »3 control the widths of the Gaussian and u' 2 controls the 
“tilt"’ of its primary axis with respect to the north direction. (By choosing this function 
we do not mean to assert that the uplift should in fact be Gaussian. We are simply using 
this method to present the observed vertical rates.) In fitting for the seven parameters 
of the model, we have used for the observation uncertainties \jo\ + of, where cq, is the 
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standard error of the uplift (see above) and cr a = 0.5 mm yr -1 . This modification is 
used to reflect the shortcomings of the simple Gaussian model. With this modification, 
the y 2 residual rate per degree of freedom was 6.5, roughly consistent with the scalings 
found in Section 3.3. 

The resulting model is shown in Figure 4b, and it can be seen that GPS-derived 
model (which we will henceforth refer to as the “observed rates”) shares much in 
common with the uplift calculated from the Earth/ice model combination (the “model 
rates”) described above. The estimated center of the uplift for the observed rates 
(A 0 = 19.5°, <p 0 = 64.2°) is quite close to, though slightly farther south than, the 
center of uplift for the model rates. The values of the observed and model maximum 
uplift rates are nearly the same, however. The areas undergoing subsidence for the 
observed field differ slightly from those for the model rates, but these areas are outside 
of the network and good agreement is not to be expected. Finally, the orientation and 
the amount of elongation agree quite well, although the model deformation is not so 
symmetric as the two-dimensional Gaussian used to represent the observed rates. This 
comparison is strong evidence that the vertical crustal motions observed with GPS are 
associated with the GIA process. 

A direct comparison between observed rates and those predicted for the Earth/ice 
model combination described above is shown in Figure 5. The excellent correlation for 
the radial velocities is clearly evident, as is the correlation for the north components. 
The east components display less of a correlation, but whether this poor agreement is 
due to greater scatter in the observed rates or an error in the ice or Earth models is not 
obvious. The observed rates show a much greater range of values, —2 to +3 mm yr -1 , 
whereas the predicted rates range only from —1 to 0.5 mm yr -1 . 

For the radial velocities, the reduced y 2 difference between the observed and model 
velocities is 14.2, or 11.6 if a mean difference of —1.0 mm yr -1 is removed. For the “n ull ” 
model that predicts zero radial velocity, the reduced y 2 is 74.3. Thus, the reduction 
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compared to the null model is a factor of 5-6. 

As a further check on our vertical rates, we compare observed sea-level rates from 
Baltic tide-gauge data to rates calculated from the Gaussian model. The tide-gauge 
data consisted of annual averages obtained from the Permanent Service for Mean Sea 
Level (PSMSL) [Pugh et al., 1987] for tide gauges with timespans of 40 years or longer 
after 1930. The exception is the tide gauge at Visby, the data for which are not in the 
PSMSL data base. The sea-level rate s at a tide-gauge located within the Baltic at (A, 
<b) is related to the land uplift u by 

s(A, <j>) = —u(X, 4>) + g(X, 4>) -I- fi (2) 

where g is the rate of change of geoid and /I is the eustatic sea-level rate. In Figure 6 
we have plotted s from tide-gauge rates versus it from our Gaussian fit. to the GPS 
vertical crustal rates. A strong correlation is evident in this figure. In effect, we have 
used our simple Gaussian model to interpolate the GPS observations to the latitude 
and longitude of the tide gauges. This correlation is clear evidence of the relationship 
between the large apparent sea-level rates observed for several centuries in the Baltic 
and the observed ongoing vertical crustal motions determined from the GPS data. These 
apparent sea-level rates have long been interpreted as indications of GIA and have even 
been used to refine Earth and ice models [e.g., Davis et al., 1999; Lambeck et al . , 1998b]. 

We conclude that the secular vertical crustal rates that we observe using the 
BIFROST GPS data are mainly associated with the ongoing GIA process. Below, we 
consider some other processes that may also contribute to the observed rates. In a 
future paper, we will present in greater detail a geophysical analysis of the BIFROST 
observations. 
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4.2. Ocean Tide Loading 

The effects of global ocean tide loading as well as the solid Earth tides are treated 
at the stage of GPS carrier phase data analysis. These motions are predominantly 
diurnal and semidiurnal; aiming for one site position estimate per day it appears more 
advantageous to account for rapid station position variations in the early processing 
stages rather than having to remove the effects a posteriori. 

The ocean loading coefficients were computed with the same method as the one 
used for the IERS Conventions 1996 [McCarthy , 1996; Schemeck , 1991]. The ocean tide 
model adopted for the processing is taken from Le Provost et al. [1994]. This model does 
not contain the Baltic Sea. a sea area that is central to our region, and the possible 
loading effects of which need a careful account and discussion. It is well established 
that the diurnal and semidiurnal tides in the Baltic Sea are less than 20 mm almost 
everywhere, and therefore their loading effects are to be expected at only sub-millimeter 
amplitudes; they can be neglected. The largest seiche-mode of the Baltic Sea. an 
east-west oscillation involving the Bay of Kiel, the Baltic Proper, and the Gulf of 
Finland occurs at 36 hr period [Wiibber and Krauss , 1979]. They can be excited by fast 
passing low-pressure areas. Significant amplitudes are found only in the bays at either 
end, lasting a couple of days. The existence of these modes argues for including either 
time-series of water level at — to the least — diurnal if not more rapid sampling rates of 
near-by tide gauges, or predicted loading effects based on such observations by means of 
a hydrodynamic model interfaced with an elastic deformation model [ Schemeck , 1991]. 

On the timescale of days to years, the situation for Baltic Sea is radically different . 
The geometry of the Baltic Sea basin being well enclosed and connected to the open 
ocean only through narrows in Denmark and between Denmark and Sweden, causes the 
mass exchange with the world ocean to be retarded. Seasonal variations of the water 
level can reach ±0.5 m as a combined response of the air pressure and wind situation, 
and due to the role of the Baltic Sea as a large catchment area where precipitation and 
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evaporation are highly variable on seasonal to interannual timescales. 

In order to obtain rough estimates of the impact of variations in the hydrology of 
the Baltic Sea area on ground deformation we have conducted a simple pilot study. 
First we assume that GPS monuments, the locations at which we aim to predict vertical 
crustal motion, are exactly following with the movement of a solid, homogeneous, 
elastic crust, i.e., we neglect porosity related deformations of soils and surface layers. 
We then devise a grid of 5 km mesh width covering the area of Sweden. Finland and 
most of Norway, on which we distinguish type of land and water coverage (a) land. 

(b) open ocean, (c) Baltic Sea, (d) great lakes, and (e) large hydropower reservoirs. 
The deformation is modeled using integrated point load Green’s functions in the usual 
way [e.g., Schemeck , 1991]. Going through all cases separately, assuming a unit height 
slab of w'ater in each of the land types, we thus can model admittance coefficients 
for the impact of loading due to (a) accumulated snow, rain and soil moisture, and 
(b)-(e) water level variations in each of the bodies. We leave a detailed account of these 
studies for future publications. Here, we use typical maximum values for the amplitudes 
of the loading processes in order to get an idea of the importance of the effects. The 
results are summarized in Table 3. 

Residual long-term rates in the water level when limiting the scope to five years 
can still be as large as 20 mm yr -1 . Thus, not accounting for these loading affects can 
offset the estimated GIA rates by 0.5 mm yr -1 . 

4.3. Atmospheric Loading 

In the EOF mixed regression we model a time-series of atmospheric loading 
for every station. Previous work on this problem [ vanDam and Wahr , 1993] showed 
generally low air pressure admittance at Onsala and small reductions of post-fit x 2 . 
The timescale of the variations in the pressure field is, unlike in the case of ocean tide 
loading, predominantly in the range of more than one day. Also, the presently available 
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processing software is not prepared for the input of three-dimensional time series of a 
priori site displacement information. 

We computed the atmospheric loading effect analogously to ocean tide loading, 
with the major difference being that we assume the loading effect is zero at the bottom 
of the open ocean assuming an inverse barometric response. Global air pressure fields 
at 1° x 1° spatial and 6 hr temporal resolution are obtained from the European Center 
for Medium-term Weather Forecasting (ECMWF) and convolved with elastic loading 
Green's functions for vertical and horizontal displacement. Although it can be shown 
that loading beyond 2000 km distance contributes little to station displacement, we 
use the entire global field. In this case the global mass balance is easy to maintain, 
and annual oscillations between the hemispheres do not offset the displacements. We 
compute an average pressure field for the entire time span in order to subtract the 
displacement due to the average atmosphere. Thus, for most of the stations a near zero 
mean for the computed pressure loading time-series is obtained. 

In the eigenvector analysis the air pressure loading information that is orthogonal 
to the station residuals is retained, while the common mode will preserve coherent 
(correlated) residual signal power from this source. Such information, however, is 
expected to be greatly suppressed since the station residuals result from a regression 
that already includes air pressure loading. 

The reason why we estimate an admittance parameter of the predicted loading 
effect rather than applying the effect as a correction is as follows. The admittance 
coefficients that we obtain are systematically and significantly lower than unity. We 
suspect that the GPS orbits induce regional perturbations since atmospheric loading is 
not applied at the stage of orbit computation. Since only a small number of stations in 
northern Europe are used by the orbit centers, only a certain fraction of displacement is 
conveyed into the orbit. 

A second issue is the possible mapping of air pressure related information into other 
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parameters of the GPS analysis. Here, the most probable candidate is the atmospheric 
delay parameter, and a particular reason to suspect it is the way the hydrostatic and 
the water vapor related delays are parameterized [Segall and Davis, 1996]. However, we 
can show that only submillimeter vertical site offsets can be expected when atmospheric 
pressure varies as much as ±30 hPa. 

Bottom pressure equilibrium requires days to weeks to establish in shallow waters. 
Since the Baltic Sea is a nearly enclosed basin, the inverse barometer response is 
certainly substantially delayed. Thus, it appears more promising to neglect the water 
response in the atmospheric loading model and add another signal channel in the linear 
regression representing the water level of the Baltic Sea at a tide gauge station nearby. 

4.4. Sea-Level Tilt 

Data from the TOPEX satellite reveals a north-south sea-level tilt during the 
duration of the TOPEX experiment (1992-present). Since this timespan is nearly the 
same as that of the BIFROST data analyzed for this study, we might expect a secular 
site position variation as a result of the elastic loading associated with this tilt. In 
order to calculate this effect, we first used the TOPEX results from the Baltic to fit 
for north-south and east-west components of the sea-level tilt. (These results were 
provided by S. Nerem of the University of Texas at Austin.) Using the Kattegatt 
to define the boundary of the Baltic, we found that within the Baltic the observed 
TOPEX sea-level rates r(A, <p) at north latitude A and east longitude q> were well 
described by r(\.<p) = A ± B { A — A 0 ) + C(o — (b 0 ), where A 0 = 18.237°, <t> 0 = 57.979°, 
A = 7.7 ± 0.2 mm yr -1 , B = —0.07 ± 0.06 mm yr -1 per degree of longitude, and 
C = —1.4 ±0.1 mm yr -1 per degree of latitude. The uncertainties were determined by 
scaling the uncertainties determined using a unit-weighted least squares solution by the 
root-mean-square (RMS) residual rate of 4.8 mm yr -1 . The distribution of the residuals 
was well described by a Gaussian distribution with the standard deviation equal to 
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the RMS residual; we did not, however, examine the residuals carefully for systematic 
geographic variations. 

In the next step, we used the tilt results to obtain a Baltic-shaped load rate with 
a north-south variation of —1.4 mm yr -1 per degree of latitude and no east- west tilt. 
We convolved this load with an elastic Green’s function in a manner similar to that 
of the other loading calculations above to calculate vertical rates in 1° steps along a 
north-south profile along longitude 20°. As expected, the maximum vertical deformation 
occurred near the south of the Baltic, at latitude 56°, where we calculated a subsidence 
rate of ~0.25 mm yr -1 . 

Thus, although the sea-level “tilt” observed by TOPEX is quite large, the overall 
load produced by the Baltic Sea, which is rather small, is nearly negligible for our 
purposes. Furthermore, one would expect that such a tilt, if caused perhaps by wind 
stress, is quite variable over longer times, so that the effective secular elastic loading 
associated with this effect will on average decrease with time. 

5. Error Analysis 

In this section we carefully examine the possible influence of a number of errors on 
our main geophysical observable, the site velocities. We include this study for several 
reasons. The expected magnitude of GIA contribution to the velocities is fairly small, 
typically < 3 mm yr -1 for the horizontal component and < 10 mm yr' 1 for the vertical. 
The technique of continuous GPS is rather new, and no careful analysis of errors has 
yet. been performed for determinations of velocity obtained from these data. We begin 
by assessing the spatial and temporal dependence of variations observed in our position 
determinations . 
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5.1. Correlation Analysis 

We anticipate that there may be a number of noise sources affecting the stations 
in a region in a similar way, such as satellite orbit errors, reference frame errors, errors 
at one or possibly more sites that are propagated through the network due to the type 
of network solution we perform, environmental conditions that change over a region in 
a coherent way (for example, soil surface reflectivity affected by climatic factors, snow 
covering antenna and radomes), and short-lived non-secular crustal deformations due to 
predominantly atmosphere and hydrosphere loading. 

Since we have time series of perturbed site positions, we seek a method that takes 
advantage of the statistics inherent in the large amount of information and relate to the 
separation of local and regional signals and noises. We have chosen to represent, the 
degree of correlation of perturbations to site position (eventually including coherent, 
transient motion) between stations using an Empirical Orthogonal Function (EOF) 
type of analysis. ( Davis and Elgered [1998] used an EOF method with estimates of 
water vapor determined from BIFROST data.) We want to utilize this information in 
the adjustment when we solve for rates, offsets and other locally relevant parameters, 
thus attempting to discriminate between local deterministic processes and correlated 
transient signals. 

Secular GIA motion is of course correlated between the sites. Thus we need a first 
stage where parameters for this process (and others we are aware of) are estimated and 
residuals formed, followed by a second stage where we then use these residuals to obtain 
improved values for the rate and other parameters. 

In the locally relevant parameter set we include atmospheric loading since this 
perturbation can be predicted on a per-site basis. However, a certain fraction of 
atmospheric loading perturbations might actually be transferred into satellite orbit 
perturbations, since corrections of site positions due to this effect are not routinely 
estimated in the precise orbit generation phase. Thus, as we will see, the atmospheric 
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loading signal is attenuated in the GPS time series, although some signal apparently 
leaks in during the second (EOF) stage due to correlation between the common mode 
and the atmospheric loading time series. 

5.1.1. Procedure. This section describes the two-stage least-squares procedure 
we have developed that uses the GPS time series of station positions to determine a 
combined set of parameters that include: (1) the usual set of rate and offset parameters; 

(2) admittance parameters for geophysical signals such as pressure and precipitation 
that we might expect to be correlated with the estimated GPS time series; and 

(3) between-site correlation parameters. We begin in the first step by looking for the 
solution to the linearized system 

G ■ p — d - f- e (3) 

The unknown parameter vector p contains both the usual set of rate and offset 
parameters and the admittance parameters for geophysical signals for one station. The 
vector d is a to-be-modeled time series (for example, north, east, or up components of 
position for one station) and e its vector of errors. The matrix G is the design matrix 
defined in the usual way. If we are including surface pressure in the model, then one of 
the columns of G will be the surface pressure at the site corresponding to the epochs for 
which the determinations of d were obtained. 

As usual, the solution to (3) is obtained by minimizing the mean square residual 
and can usually be determined using the standard least-squares formulation. However, 
since this system can in general be underdetermined or otherwise singular, and for 
computational compatibility with the second step of the process, we solve (3) using the 
generalized inverse method of Lancxos [Aki and Richards, 1980]. In this method, we first 
normalize the system of equations (3) such that the expected variance of the e is the 
unit matrix. (We assume for this analysis that the observations are uncorrelated.) We 
denote the normalized quantities with a tilde (G, d, and so on). The generalized inverse 
of G is expressed as a singular value decomposition (SVD), wherein we seek the solution 


29 


to the eigenvalue equations 



where it is a matrix of eigenvectors that span the parameter space, v is a matrix of 
eigenvectors that span the data space, and A represents the set of eigenvalues. Aki and 
Richards [1980] demonstrate that u and v share the first m eigenvalues, where m is the 
minimum of the dimensions of u and v\ the remaining eigenvalues are zero. 

In our analysis for this study, w r e modeled the time series of vertical components of 
station position for the BIFROST sites. In addition to rates and biases for each site, 
we estimated site-dependent sinusoidal amplitudes (in- and out-of-phase) with annual, 
semiannual, terannual, and quarterannual periods, and (site-dependent) admittance 
parameters for surface pressure. 

In the second stage of the analysis the solution was repeated with the modification 
that, for each site we also estimated parameters that represent admittance for the 
residuals calculated in the first stage for all the BIFROST sites. For example, if using 
the data from the Umea site yields a non-zero admittance for the stage-one residuals 
from Kiruna site, that result indicates a correlation between the unmodeled signal 
from the first stage for Kiruna and the observed signal for Umea. In this sense, the 
second stage is equivalent to a traditional (unweighted) EOF analysis. Moreover, in a 
manner analogous to a traditional EOF analysis, we keep only those parameters wdiose 
eigenvalues signify that they convey a significant amount of information. 

In the second (EOF) stage, we define two subsets of the parameter space. The 
subset of the parameters that were estimated in the first stage we name C: these 
parameters will be estimated regardless of their eigenvalue. The subset consisting of the 
new' (residual admittance) parameters of the EOF construct we name E. Sorting the 
eigenvalues of £ by decreasing magnitude yields one w'hich is the largest, A c . Associated 
with it is eigenvector v c , which is to be retained; the remaining subspace of £ is to be 
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ignored. In a traditional space-time EOF analysis v c is known as the first temporal 
eigenvector. 

The eigenvector v c retained from subspace £ extends the data space to comprise 
a common mode. The common mode time-series is obtained from the parameter 
solution by simply projecting it on the subspace £ and inserting it into the model (3). 
The common mode construction is a form of a spatial filter performing a weighted 
mean, where the weights are devised by the eigenvalue process and hence by the data 
themselves. The procedure we devised may be contrasted to the filter of Wdowinski et 
al. [1997], which was formed by taking a mean of the time series for the stations in the 
network. Each observation is equally weighted regardless of uncertainty or station. In 
our procedure, we use the uncertainties and the station weights are determined by the 
station admittance parameters, calculated to minimize the mean-square residual. 

As in a traditional EOF analysis, the relative signal power propagated through 
each eigenvector in the inverse solution is proportional to the squared magnitude of 
the eigenvalue. In our solutions using the vertical component of site position, the 
ratio X 2 / A 2 is typically on the order of 20%. The normalized x 2 residual for the 
second stage is typically a factor of ~2 smaller than that of the first stage, wherein no 
common-mode signal is estimated. In Figure 7 we show the normalized residuals for the 
first and second analysis stages along with the common mode signal. 

5.1.2. Spatial correlations. The eigenvectors and eigenvalues carry information 
about the correlation of the time series among the sites, accounting for the model of 
the first adjustment stage. This information stems from the similarity of the scalar 
product used in the eigenvalue solution process and the manner in which correlations 
are computed. Figure 8, for example, contains a plot of the parameter eigenvectors for 
an analysis of a subset of 11 SWEPOS sites: Metsahovi, Skelleftea, Umea, Sundsvall, 
Martsbo, Leksand, Sveg, Ostersund, Vilhelmina, Arjeplog, and Kiruna. The 11 
admittance parameters representing £ are the last 11 (parameter numbers 15-25). It is 
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clear from the primary eigenvector in Figure 8 that the contribution from the different 
sites is fairly large and constant. The exception might be the first site, Metsahovi. 
This site is used to constrain the reference frame (see above). The correlation between 
regional stations and a station which has been used for mapping into the reference frame 
is generally low because the average rotations and translations of the the small number 
of constraining stations that have been used are only slightly overdetermined; thus they 
are propagated throughout the network. 

It is straightforward to show that the correlation coefficient 7 j* between time series 
from the jth and kth stations can be computed using the first spatial eigenvector of £. 
denoted u c ; 

(h) 

where N is the number of data. The expected low correlations with constraining 
stations, in this caseMetsahovi (parameter #15). is clearly seen in Figure 9, which 
shows the correlation as a function of intersite distance. The remaining rate parameters 
exhibit great coherence (7 ~ 0.5) on a regional scale. Figure 9 does seem to indicate, 
however, that there is a clear but weak dependence of correlation on intersite separation. 
This result implies that the cause of the correlation is network-wide, indicating perhaps 
a reference frame or orbital- type effect. 

5.1.3. Effect of reference frame errors. There are two types of reference frame 
errors to consider: errors formally accounted for by the error covariance matrix, and 
biases. In the ideal case of independent Gaussian measurement errors with perfect 
geodetic models, linear propagation of the measurement error covariance would yield 
an accurate characterization of the errors in geodetic model parameter estimates. 
Uncertainties in velocity estimates derived from such position determinations would 
similarly be accurate and well understood. In the real world, however, geodetic 
models are not perfect, measurement noise processes are not normally (i.e., Gaussian) 
distributed, and the measurement errors are not generally independent. Each of 
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these factors complicates the parameter estimation problem. The resulting parameter 
estimates are likely to be contaminated by systematic errors which are difficult, perhaps 
impossible, to fully assess. 

An exhaustive list of the mechanisms through which systematic errors could 
manifest themselves as reference frame biases and their respective importance will 
require continued research. However, it is not difficult to list some important potential 
mechanisms. Deficiencies in the geodetic models describing satellite orbital dynamics 
or the dynamics of Earth's orientation, for example, have obvious implications for an 
accurate reference frame realization. For local networks, errors in the atmospheric 
modeling or other spatially correlated errors, such as those due to similar scattering at 
similar monuments, could result in local reference frame biases. The same is also true of 
over-constrained a priori parameter estimates, such as for satellite ephemerides, satellite 
clock variations, or site positions, particularly if these a priori estimates were themselves 
correlated. These “fiducial errors” have been appreciated for some time [e.g., Larson 
and Agnew. 1991]. Site specific errors, associated with such phenomena as multipath or 
antenna phase center variations, may manifest themselves indirectly as reference frame 
biases as we now discuss. 

Trade-offs between parameters of the geodetic model, such as between the implicit 
specifications of the satellite and GPS network orientations, render the site position 
determination problem ill-posed. That is, in the absence of additional constraints, the 
matrix of partial derivatives relating the basic prefit GPS carrier phase and pseudo range 
residuals (i.e., the data) to the first order corrections to the a priori parameter estimates 
(including site positions, satellite parameters, and Earth orientation parameters) is 
singular. The degree of singularity depends on the geometry of the tracking network; 
regional scale networks are “less singular” than global scale networks. A priori 
information, such as knowledge of the positions of certain sites in the network, is often 
effectively employed to regularize this singularity. The resulting position estimates 
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implicitly define a reference frame. Depending on the level of constraints imposed and 
assuming, for now, that this level is consistent with the true errors in the a priori 
estimates, some linear combinations of the model parameter estimates, such as the 
differences in positions of the sites, may be significantly better determined than the 
‘'absolute” model parameters themselves. However, site specific errors not formally 
accounted for by the covariance matrix will affect all components of the vector of linear 
combinations involving the contaminated site. The importance of any reference frame 
errors resulting from this transformation will thus depend on the transformation itself 
(the particular linear combination), the number of sites suffering from systematic errors, 
and the size of these systematic errors. 

5.2. Power Spectral Analysis 

The time series contains signals that are not well modeled by simple rate and 
offset terms. To assess these errors it is most useful to examine their power spectrum. 
Generally, one computes a set of postfit residuals, then uses the residuals as input to. for 
instance, an FFT algorithm. Windowing or filtering may be performed at some point. 
In almost all cases the algorithm for computing the power spectrum assumes that the 
data are equally spaced and equally weighted. Given the nature of our data, neither of 
these assumptions is very good. Moreover, it would be useful to be able to assess the 
effects of data gaps and the estimation of model parameters, especially rate and offset 
parameters, on the estimated power spectrum. 

We have implemented a method for estimating a limited set of components for 
the power spectrum of the unmodeled signals in our time series. In this method, we 
estimate a series of coefficients to sinusoidal terms simultaneously with the rate and 
offset parameters. In effect, we used a very truncated form of this when we estimated 
coefficients to annual sinusoids (see above). The following discussion will apply to data 
from sites with long time series only, the SWEPOS sites and Metsahovi. 
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Given the total timespan of four years for these series of daily estimates, an 
FFT-tvpe algorithm would estimate spectral components for frequencies 
k = where A/ min = 0.25 cycles per year (cpy) and N ~ 1460. It is the 

lowest- frequency components that are of most interest, since the presence of power 
at these frequencies might indicate errors having correlation times long enough to 
significantly influence our rate estimates. However, because we are estimating a rate 
parameter and offset parameters that can absorb power at these lower frequencies, we 
might expect that solutions containing low-frequency power-spectrum terms might be 
fairly unstable. In fact, we see exactly this result. Although it depends on the specific 
site, in general we cannot resolve well the power-spectrum parameters for frequencies 
lower than about 1 cpy. 

We have therefore included spectral components for frequencies kAf min with 
A fmi-n — 1 cpy and 1 < k < 64. Above this maximum frequency, the spectral 
components are generally not significant. Examples of the power spectra are shown in 
Figure 10. The spectra tend to be reddened, especially those for the vertical component 
of site position. The spectra for the vertical components are greater in magnitude 
than those for the horizontal components, as might be expected. The spectra for the 
horizontal components tend to be flatter. The north component especially tends to have 
(theoretically) significant peaks at higher frequencies. 

The spectra we have calculated indicate that there is a great deal of low-frequency 
power in the variations of site position that is not properly modeled as a secular 
variation with offsets. The source of these variations cannot be determined from the 
power spectra alone. Spectra of residuals from our EOF analyses (not shown) indicate 
that some of the power is associated with network-wide effects. Some, but not all, of the 
power is therefore associated with site or region-wide sources of error. For example, we 
know that the snow effects (see below) cause errors with pseudo- annual periods. From 
our analysis of monument stability (see below), it seems unlikely that a large component 
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of the power is associated with highly localized monument motions. 

Visual inspection of the power spectra lead us to conclude that the spectral indices 
for these spectra are in the range of —1 to —2, consistent with the analyses of Zhang 
et al. [1997] for GPS sites in southern California and Mao et al. [1999] for global sites. 
Again, the horizontal components tend to be a little flatter, with spectral indices closer 
to —1. As pointed out by these authors and others, if these spectral components are 
representative of a stochastic noise process with a long correlation time, then the 
standard errors one determines for the slopes are significantly underestimated. However, 
it is probable, as we have pointed out, that our error spectra contain some truly 
periodic terms. Indeed, many of the spectra shown by Mao et al. [1999] have peaks 
near annual frequencies. It is clear, though, from our attempts to estimate spectral 
components with frequencies lower than 1 cpy, that noise at these lower frequencies 
would significantly effect our estimates of the rates. Further studies regarding the 
detailed spectral composition of our time series will be possible when we have longer 
timespans of determinations with no changes in equipment. 

5.3. Atmosphere 

As discussed in Section 3, the GIPSY analysis software uses a stochastic filter 
to model the temporal variability of the zenith atmospheric propagation delay. 

The a priori model for the propagation delay for each site consists of a constant 
ellipsoidal-height-dependent term representing the hydrostatic zenith delay [Davis et al.. 
1985] that is "mapped” to the elevation angle of the GPS signal using the Lanyi dry 
mapping function [Lanyi. 1984]. The a priori wet zenith delay is taken to be 100 mm. 
Corrections to this a priori total zenith delay are estimated by adopting that the Lanyi 
wet mapping function to describe how the correction maps with elevation angle. In 
effect, this procedure maps the combined corrections to both the hydrostatic and the 
wet zenith delays using the wet mapping function. Using a simulation procedure similar 
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to that of Eldsegui et al. [1995], we find that the error in the vertical component of site 
position (the primary geodetic parameter affected) caused by mapping the correction to 
the hydrostatic zenith delay using a wet mapping function is less than ±2 mm, assuming 
an annual pressure variation of ±30 mbar. 

The mapping functions used assume that, the atmosphere is azimuthally symmetric 
as viewed from the station. Horizontal structure in the atmosphere nevertheless 
generally leads to an effective azimuth-dependent, mapping function. In particular, 
horizontal gradients in the refractive index lead to sinusoidal variations in azimuth 
[e.g., Davis et al., 1993]. We have not incorporated this model into our phase model or 
estimated any parameters associated with horizontal structure, leaving this to future 
studies. Although such errors can be quite large at low elevation angles [Davis et al., 
1993], we do not expect this error to bias systematically the estimate of the site rates. 
Of much greater concern are the azimuthal-dependent effects of snow accumulation on 
the radomes. 

5.4. Signal Effects 

In the course of our experiment we have investigated several errors involving 
the propagation of the GPS signals in the nearby environment of the GPS antenna. 
Near-field scattering [Eldsegui et al., 1995] involves the reflection of the GPS signal off 
surfaces close to the GPS antenna. This effect is different from multipath in that the 
reflecting surfaces are well within the near-field of the antenna. These surfaces thus 
become electromagnetically coupled to the antenna, effectively creating a new antenna. 
One such reflecting surface is the pillar directly beneath the GPS antenna, which often 
(and in the case of our Swedish sites) has a metal plate embedded onto its top surface. 
Jaldehag et al. [1996b] investigated this effect for sites of our Swedish GPS network, 
and found that for the most part the scattering effects are equal for all the sites (and 
therefore cancel in the baselines), since with the exception of the Onsala GPS site all 
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the pillars are identical and the local vertical direction at each of the sites is nearly 
parallel. It is possible that this effect is in part responsible for the common-mode site 
variations we find (Section 5.1). In the principal component analysis we performed, we 
found that the component of the spatial eigenvector for the Onsala site had the smallest 
value of all the Swedish sites. 

The error associated with the signal scattering effect is quite large, even when the 
size of the phase error is small. In the least-squares estimation procedure, the errors 
in each of these parameters can be “magnified" and yet the summed contribution to 
the phase can be fairly small. The effect on the horizontal components is smaller than 
that for the vertical but not zero. Eldsegui et al. [1995] found that the error in the 
estimated vertical component was therefore dependent on the geometrical distribution 
of observations on the sky. For a given fixed constellation and sequence of observations, 
the error in the estimate of site position is dependent on the minimum elevation angle 
accepted for use in the data analysis. In practice, however, it is not possible to maintain 
a fixed geometrical distribution of observations on the sky, due to data dropouts. 

We found also that snow and ice, which adheres to the radome and accumulates, is a 
significant source of error, simply by causing an additional propagation delay. This effect 
was first noted by Webb at al. [1995] when several meters of snow buried a GPS antenna 
near Long Valley. The amount of snow and ice which accumulates on our antennas is 
much smaller, perhaps several decimeters maximum, and the problem is therefore more 
insidious. (The lower amount of snow also tends to be distributed unevenly on the 
radome.) Jaldehag et al. [1996a] found that they were able to approximate the results 
from elevation angle cutoff tests using a conically symmetric distribution of snow, but 
visual inspection of our sites indicates that real conditions are not so ideal. 

A final effect which we have already mentioned is the propagation delay due to 
the radome. In theory an expression can be developed for this contribution and in the 
future we plan to test such a model. For this paper, however, we have simply introduced 
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an offset to the time series on the epochs which we have changes radomes. 

5.5. Monument Stability 

As described in the section describing the GPS networks, we monitor the relative 
position of the pillar reference point within a local (10-15 m area) network of reference 
pins. The determination of the pillar position is given in Table 4. The expected 
uncertainty of these determinations, based on propagation of the instrumental errors 
through the resection technique, is ~0.1 mm. Based on the RMS analysis of all the 
data, the standard deviation of a single measurement is 0.2 mm for the north, 0.3 mm 
for the east, and 0.4 mm for the vertical. However, there are several occasions when 
repeated measurements are obtained in a single day. The five repeated measurements 
from the Leksand site, for instance, yield RMS differences of 0.01 mm for north, 0.04 for 
east, and 0.06 mm for vertical. The single repeated measurement for Kiruna. on the 
other hand, yields an RMS difference (based on two measurements) of 0.2 (north), 

0.3 (east), and 0.1 mm (up), roughly consistent with the same statistics for the single 
repeated measurement for the Overkalix site: 0.03 (north), 0.1 (east), and 0.1 mm (up); 
and for the Norrkoping site: 0.1 (north), 0.1 (east), and 0.1 mm (up). If we combine 
the Kiruna, Overkalix, and Norrkoping repeat measurements, we find an RMS error 
of 0.10 (north), 0.18 (east), and 0.10 (north). We will therefore adopt the average 
0.13 mm as our uncertainty in each of the components. We have no explanation as to 
w r hv the differences in the repeated Leksand measurements are so small relative to these 
values. Eliminating the multiple measurements by replacing them with their average, 
the overall RMS variations for all sites are 0.20 (north), 0.31 (east), and 0.38 mm (up). 
Thus, overall, the average RMS variations axe consistent with the adopted measurement 
error at the level of 3 a. 

Nevertheless, for individual pillars, most notably Skelleftea S (east component) 
and Umea (up component), the variations are greater than 1 mm. In the case of 



39 


Skelleftea S, the variation for the north and vertical component are small, so unless the 
pillar happened to have moved on a nearly north-south axis this result is spurious. The 
Umea variation is in the vertical component, which, one can surmise, could experience 
motions independent from the horizontal components. The most worrisome effect would 
be settling of the pillar. However, the measurements indicate that the pillar was 2.3 mm 
higher relative to the nearby network in June 1995 than in June 1993. In order to resolve 
this and other issues, further pillar measurements are planned for the near future. 
However, these future measurements will be obtained in a slightly different manner. 
Reflectors will be fixed to the pillars and the antenna will not have to be removed. 

More frequent measurements were obtained on the Leksand site than for others, 
for test purposes. The time series of measurements is shown in Figure 11. Langbein 
and Johnson [1997] have investigated temporal variations in line lengths determined 
from EDM; their analysis indicates that the PSD of the variations follows a f~ 2 
behavior, indicative of a random walk process. They argue that these variations are 
associated with pillar motions. Although our Leksand data are sparse, we can perform a 
maximum-likelihood analysis on the assumption that the temporal variations of the pillar 
position can be described by a random walk, and that the measurement uncertainties 
are 0.13 mm, as described above. For the three components, the maximum likelihood 
estimate of the random walk variance was 0.3 (north), 0.2 (east), and 0.1 mm 2 yr -1 
(up), values much less than the average value of 1.7 mm 2 yr -1 found by Langbein 
and Johnson [1997], The monuments investigated by these authors are very similar to 
those we used for the Swedish GPS network: for instance, the sites of the Parkfield 
network used galvanized steel pipe of diameter 20 mm placed into augered holes 2 m 
deep. Near two of these Parkfield sites, special monuments were installed to a depth of 
10 m; the average random walk variance for these two deep anchored monuments was 
~2 mm 2 yr -1 , whereas that for the shallow monuments was ~15 mm 2 vr _1 [Langbein 
and Johnson, 1997], 
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The average random walk variances, based on the difference between end-point 
measurements for the 20 sites with pillar measurements spanning more than six months, 
is 0.08 (north). 0.21 (east), and 0.26 mm 2 yr -1 (up). 

6. Discussion and Future Work 

We have used a regional continuous GPS network to measure three-dimensional 
crustal deformation rates in Fennoscandia. The observed rates correlate highly to rates 
predicted for a viscoelastic Earth undergoing present-day glacial isostatic adjustment 
(GIA) due to the rapid melting of Late Pleistocene glaciers. The observed rates correlate 
well also with observed sea-level rates obtained from tide-gauge data over the last 30 
years. 

We see no evidence for the regional shear inferred by Pan and Sjoberg [1999] using 
a much smaller GPS network and only two campaign measurements. Pan and Sjoberg 
[1999] found relative horizontal motions across the north-south extension of the Baltic 
Sea to be 2-3 mm yr -1 , whereas using the Table 2 rates we find the motions of these 
sites to be generally less than 1 mm yr -1 , probably within the uncertainties. Our 
velocities are roughly consistent with the expected velocities near the center of uplift . 

We get much better agreement with theoretical predictions of crustal velocity based 
on GIA, relative to the “standard errors,” in the radial (vertical) component than in the 
horizontal components of station velocity. This results may seem counterintuitive since 
radial velocity estimates obtained from GPS data are well known to be less accurate that 
horizontal velocity estimates. Nevertheless, there could be several explanations for this 
situation. The standard errors for the vertical component could be more realistic, since 
they are by nature larger. If, for example, there were some systematic error (correlated 
in time or not) that was about equal in its effects on the three components of site 
position, then proportional to the standard error the systematic error would be much 
larger in the horizontal components. Errors in the ice model would also tend to have a 
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proportionally greater effect on the horizontal components of velocity. For example, a 
slightly misplaced maximum ice thickness could easily change the predicted horizontal 
rates by 100% near the maximum, while changing the radial rate by only ~10%. 

Although we have much evidence for time-correlated behavior in the site-position 
time series both from power spectrum calculations and from visual inspection of the 
time series themselves, our measurements for pillar motion indicate that the top of the 
pillar moves less than 1 mm relative to a very local network of pins. This result may 
indicate that “monument instability/’ if it exists for our sites, may not be the result of 
local deformation of the crust. If so, then the “footprint” would have to be extended in 
order to measure this phenomenon. 

The empirical orthogonal function (EOF) analysis we performed revealed a 
correlated noise in radial position time series across a large region (~1000 km). This 
analysis also indicated that these correlations decrease slightly with distance. Typical 
RMS residuals for our radial time series are 7 mm. The EOF analysis indicated that 
over half of the noise contributing to the RMS residual is from this correlated noise. The 
source of the correlated noise, unfortunately, cannot be revealed by an EOF analysis. 
We and others have speculated that the correlations are due to systematic errors in the 
satellite orbits, terrestrial reference frame, or both. If the source of this error can be 
identified and eliminated, then there is the possibility of reducing velocity uncertainties 
from GPS by ~50%. 

In the near future, we will also extend our analysis to include more recent data. It 
is clear from our analysis that the several changes we have made to the GPS equipment 
have limited the achievable accuracy. Longer time series with no equipment changes 
should provide more reliable rate estimates. 

We have not attempted to determine our “true” uncertainties. We have reported 
the “standard errors" so that users of the rates may have the straightforward results of 
the least-squares analysis. The standard errors may be scaled or otherwise increased as 
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required. Based on a linear fit to the predictions we presented in the paper, the standard 
errors would have to be increased a factor of ~8 for the horizontal components and a 
factor of ~3 for the vertical to achieve y 2 values of unity. These are much larger than 
the scale factors found in Section 3.3. The comparison to the model assumes, however, 
that all the disagreement between the observed and predicted rates is caused by errors 
in the observed rates. In a future work, we will use the difference between the observed 
and predicted rates to infer errors in the the Earth and ice models used to obtain the 
predictions. The degree of fit can then be used to asses much better the errors in the 
observed rates. 

The results presented here provide the first dense regional geodetic determinations 
of ongoing crustal deformation associated with GIA. The rates estimated herein are 
intended to be used as a new observable in the study of this phenomenon. We will 
also use the radial rates to correct tide-gauge data for vertical crustal motion to obtain 
absolute sea-level rates. 
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Table 1 . BIFROST site histones. 

Approximate Position 0 



IERS Net- 

Installed/ 



Ra- 


North 

East 

Height. 

Site 

Domes No. 6 work c 

Modified 

Receiver^ 

Antenna e 

dome-f 

Monument 9 

Latitude 

Longitude 

m 

Arjeplog 

S 

93/08/20 

8000 

DM T 

Delft 

SWEPOS 

66° 19' 

18° 07' 

489.2 



95/02/03 

95/08/01 

Z-XII 


Type A 







96/06/28 



None 







96/10/29 



Type B 





Boras 

S 

96/11/12 

Z-XII 

DM T 

Ash 

SWEPOS 

57° 43' 

12° 53' 

220.0 

Hassleholm 

S 

93/08/20 

8000 

DM T 

Delft 

SWEPOS 

56° 06' 

13° 43' 

114.1 



95/05/19 

95/06/10 


* 

Type A 







95/08/01 

96/07/01 

Z-XII 


None 







96/11/12 



Type B 





Joensuu 

10512M001 F 

95/06/15 

Z-XII 

DM Ash 

Ash 

2.5 m SG 

62° 23' 

30° 06' 

113.7 

Jonkoping 

S 

93/08/20 

8000 

DM T 

Delft 

SWEPOS B 

57° 45' 

14° 04' 

260.4 



94/06/23 

95/06/24 


* 

Type A 







95/05/21 

95/08/01 

96/07/03 

Z-XII 

* 

None 







96/11/10 



Type B 





Karlstad 

S 

93/08/20 

8000 

DM T 

Delft 

SWEPOS 

59° 27' 

13° 30' 

114.3 



95/02/08 

95/05/23 

95/08/01 

Z-XII 

* 

Type A 







96/07/03 



None 







96/11/15 



Type B 





Kevo 

F 

96/07/05 

Z-XII 

DM Ash 

Ash 

5 m SG 

69° 45' 

27° 00' 

135.9 

Kiruna 

S 

93/08/01 

8000 

DM T 

Delft 

SWEPOS 

67° 53' 

21° 04' 

498.0 



93/08/17 


* 








94/06/15 


* 








95/06/16 

95/07/10 

95/08/01 

Z-XII 

* 

Type A 







96/07/30 



None 







96/10/30 



Type B 





Kivettv 

F 

96/03/05 

Z-XII 

DM Ash 

Ash 

2 m CP 

62° 49' 

25° 42' 

216.3 

Kuusamo 

F 

96/10/30 

Z-XII 

DM Ash 

Ash 

2.5 m SG 

65° 55' 

29° 02' 

379.0 

Leksand 

S 

93/08/01 

8000 

DM T 

Delft 

SWEPOS 

60° 43' 

14° 53' 

478.1 



93/08/11 


* 








94/01/18 


* 








94/02/06 


★ 








94/03/08 


* 








94/04/15 


* 








94/06/14 


♦ 








94/08/25 

95/01/30 

95/08/01 

Z-XII 

* 

Type A 







95/10/05 

96/06/27 


* 

None 







96/10/25 



Type B 





Lovo 

S 

93/08/01 

8000 

DM Ash 

Delft 

SWEPOS C 

59° 20' 

17° 50' 

79.7 



93/10/28 

95/05/16 


* 

Type A 
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95/06/15 
95/08/01 Z-X1I 
96/06/28 
96/11/07 

* 

None 
Type B 





Metsahovi 

10503S01 1 

F 

92/01/01 SNR-C 
95/04/30 8100 

DM B 

None 

25 m SG IS 

60° 13' 

24° 24' 

94.6 

Marts bo 


S 

93/08/01 8000 
95/02/07 
95/08/01 Z-XI1 
96/05/07 

DM T 

Delft 
Type A 

Type B 

SWEPOS C 

60° 35' 

17° 16' 

75.4 

Norrkoping 


s 

93/08/01 8000 

95/05/15 

95/08/01 Z-XII 

95/10/06 

96/07/12 

96/11/09 

DM Ash 
* 

Delft 
Type A 

None 
Type B 

SWEPOS 

58° 35' 

16° 15' 

41.0 

Olkiluoto 


F 

94/10/19 8100 
95/04/19 Z-XII 

DM T 

Delft 

2 m CP 

61° 14' 

21° 28' 

30.5 

Onsala 


S 

93/07/01 8000 

93/08/16 

95/05/20 

DM B 

Delft 

IGS 

57° 24' 

11° 56' 

45.5 

Oskarshamn 


s 

93/08/01 8000 

95/05/18 

95/06/13 

95/08/01 Z-XII 

96/06/29 

96/11/11 

DM Ash 
* 

Delft 

Type A 

None 
Type B 

SWEPOS 

57° 04' 

15° 60' 

149.8 

Oulu 


F 

95/09/16 Z-XII 

DM Ash 

Ash 

8 m IS 

65° 05' 

25° 54' 

79.5 

Romuvaara 


F 

96/05/07 Z-XII 

DM Ash 

Ash 

2 m CP 

64° 13' 

29° 56' 

241.7 

Skelieftea 


S 

93/08/01 8000 

93/08/15 

95/02/24 

95/06/15 

95/08/01 Z-XII 

96/07/03 

96/11/12 

DM T 
* 

* 

Delft 
Type A 

None 
Type B 

SWEPOS 

64° 53' 

21° 03' 

81.2 

Sodankylla 

10513M001 

F 

94/08/14 8100 
95/05/15 Z-XII 

DM T 

Delft 

2.5 m SG 

67° 25' 

26° 23' 

299.8 

Sundsvall 


S 

93/08/01 8000 

95/02/06 

95/06/13 

95/08/01 Z-XII 

96/07/01 

96/11/04 

DM T 
* 

Delft 
Type A 

None 
Type B 

SWEPOS 

62° 14' 

17° 40' 

31.8 

Sveg 


s 

93/08/01 8000 

95/01/31 

95/06/21 

95/08/01 Z-XII 

96/07/01 

96/10/26 

DM T 
* 

Delft 
Type A 

None 
Type B 

SWEPOS 

62° 01' 

14° 42' 

491.2 

Tuorla 


F 

94/08/15 8100 
95/01/21 Z-XII 

DM T 

None 

2.5 m SG 

60° 25' 

22° 27' 

60.5 

Umea 


S 

93/08/01 8000 
95/02/05 
95/06/14 
95/08/01 Z-XII 

DM T 
* 

Delft 
Type A 

SWEPOS 

63° 35' 

19° 31' 

54.5 
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96/08/13 

96/11/03 



Vaasa 10511M001 

F 

95/04/07 

Z-XII 

DM Ash 

Vilhelmina 

S 

93/08/01 

95/02/02 

8000 

DM T 



95/06/18 

95/08/01 

96/06/28 

96/10/27 

Z-XII 

* 

Virolahti 

F 

94/08/15 

8100 

DM T 



95/03/24 

Z-XII 


Vis by 

S 

93/08/01 

95/06/14 

8000 

DM T 



95/08/01 

96/06/25 

96/11/08 

Z-XII 


Vanersborg 

s 

93/08/01 

8000 

DM T 



93/09/09 


* 



95/05/22 

95/06/24 

95/08/01 

96/05/23 

96/11/13 

Z-XII 

* 

Ostersund 

s 

93/08/01 

95/02/01 

8000 

DM T 



95/08/01 

95/09/15 

96/07/08 

96/10/27 

Z-XII 

* 

Overkalix 

s 

93/08/01 

8000 

DM Ash 



94/06/15 


* 



95/06/16 

95/07/10 

95/08/01 

96/06/28 

96/11/01 

Z-XII 

* 


None 
Type B 





Ash 

2.5 m SG 

62° 58' 

21° 46' 

58.0 

Delft 
Type A 

SWEPOS 

64° 42' 

16° 34' 

450.0 

None 
Type B 





Delft 

2.5 m SG 

60° 32' 

27° 33' 

36.9 

Delft 
Type A 

SWEPOS 

57° 39' 

18° 22' 

79.8 

None 
Type B 





Delft 

SWEPOS 

58° 42' 

12° 02' 

169.7 

Type A 





None 
Type B 





Delft 
Type A 

SWEPOS 

63° 27' 

14° 51' 

490.1 

None 
Type B 





Delft 

SWEPOS 

66° 19' 

22° 46' 

223.0 


a WGS-84 ellipsoidal coordinates. 

b The DOMES number is a unique station identifier issued by the International Earth Rotation Service. Only some 
of the SWEPOS sites have been so registered. 

C S = SWEPOS, F = FinnRef. 

d 8000 = AOA SNR-8000; 8100 = AOA SNR-8100; Z-XII = Ashtec Z-XII. 

e Variants of the Dorne-Margolin (DM) chokering GPS antenna have been used, including the DM-B, the DM-T, 
and the Ashtec manufactured version denoted ‘DM Ash.” An asterisk indicates that the antenna was removed and 


Type A 

None 
Type B 


replaced. 

f See text and Emardson et al. [2000]. 

9 SWEPOS Monuments are denoted as “SWEPOS” for standard, and variants “B” and “C” (see text). IGS is the 
Onsala mount (see text). SG = Steel Grid Mast. CP = Concrete Pillar, IS = Invar Stabilized. 
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Table 2. Site velocities in ITRF96. 


Velocity Estimates, mm yr 1 

Standard Solution Edited Solution EOF Solution 


Site 

East 

North 

Up 

East 

North 

Up 

Up 

Arjeplog 

— 1.5±0.1 

1.4±0.1 

10.2±0.8 

— 2.0±0.3 

1.4±0.2 

10.1±0.8 

9.2±2.3 

Boras 

— 2.0±0.2 

-O.liO.l 

1.5± 1.1 

— 2.0±0.5 

— 0.1±0.3 

2.9± 1.3 


Hassleholm 

1.0±0.1 

0.3±0.1 

2.6±0.6 

— 0.3±0.2 

— 0.8±0.1 

— 1.5±0.6 

3.4± 1.3 

Joensuu 

1.7±0.2 

-0.6±0.1 

— 1.5±0.4 

1.1±0.2 

— 0.4±0.1 

— 1.1±0.5 


Jonkoping 

l.liO.l 

-0.2±0.1 

0.7±0.8 

— 2.3±0.3 

— 0.8±0.2 

0.2±0.9 

1.7± 1.7 

Karlstad 

— 0.5±0.1 

-0.5±0.1 

4.3±0.6 

— 1.4±0.3 

— 0.4±0.1 

2.9±0.7 

4.1± 1.3 

Kevo 

1.8±0.3 

1.7±0.2 

— 3.5±1.1 

2.1±0.3 

1.5±0.2 

— 0.7± 1.1 


Kiruna 

0.3±0.1 

1.5±0.1 

11.7±0.9 

0.3±0.3 

— 0.3±0.2 

4.7±0.9 

10.6± 1.9 

Kivetty 

— 1.3±0.4 

— 0.8±0.2 

5.5± 1.0 

0.2±0.4 

— 0.8±0.2 

3.2± 1.0 


Kuusamo 

2.9±0.5 

— 0.2±0.3 

8.8±1.4 

2.9±0.5 

— 0.4±0.3 

2.9±1.6 


Leksand 

-O.liO.l 

0.7±0.1 

5.1±0.7 

— 1.5±0.4 

0.2±0.2 

4.6±1.0 

5.2±1.8 

Lovo 

1.5±0.1 

— 0.7±0.1 

1.5±0.9 

-0.8±0.4 

— 0.4±0.2 

0.6±1.0 


Metsahovi 

O.liO.l 

— 1.5±0.1 

2.1±0.7 

— 0.3±0.1 

— 1.6±0.0 

1.4±0.2 

3.7± 1.4 

Martsbo 

— 0.7±0.1 

— 0.3±0.1 

6.5±0.6 

— 0.5±0.2 

-1.0±0.1 

6.0±0.6 

5.7±1.4 

Norrkoping 

1.7±0.1 

-0.8±0.1 

3.2±0.9 

— 0.6±0.4 

-0.8±0.2 

1.9± 1.0 


Olkiluoto 

1.9±0.2 

-1.0±0.1 

7.5±0.5 

1.4±0.2 

— 0.6±0.1 

7.0±0.5 


Onsala 

— 0.2±0.1 

— 0.4±0.1 

0.1±0.8 

— 1.4±0.3 

— 0.4±0.1 

— 1.0±0.7 

0.3±1.5 

Oskarshamn 

l.liO.l 

0.8±0.1 

1.1±0.9 

— 1.6±0.4 

— 0.3±0.2 

-0.8± 1.0 


Oulu 

0.7±0.2 

0.2±0.1 

5.9±0.6 

1.2±0.2 

0.2±0.1 

3.7±0.6 


Romuvaara 

— 0.2±0.5 

— 1.5±0.3 

7.7± 1.5 

— 1.8±0.6 

— 0.9±0.3 

11.4± 1.6 


Skelleftea 

— 0.3±0.1 

O.liO.l 

9.8±0.8 

0.2±0.3 

— 0.1±0.2 

8.1±0.9 

10.7± 1.6 

Sodankyla 

1.2±0.2 

0.6±0.1 

6.1±0.5 

1.2±0.2 

0.6±0.1 

6.5±0.6 


Sundsvall 

— 1.4±0.1 

O.liO.l 

8.6±0.6 

— 1.9±0.3 

— 0.3±0.1 

7.2±0.7 

8.8± 1.6 

Sveg 

— 0.4±0.1 

0.2±0.1 

7.5±0.6 

— 1.7±0.3 

-O.liO.l 

7.7±0.7 

9.9± 1.6 

Tromso 

— 1.8±0.1 

1.5±0.1 

— 0.8±0.3 

— 1.5±0.1 

1.7±0.1 

— 0.7±0.3 


Tuorla 

1.3±0.2 

-0.6±0.1 

2.7±0.5 

1.3±0.2 

— 0.6±0.1 

3.4±0.5 


Umea 

— 0.2±0.1 

O.liO.l 

10.1±0.6 

— 1.1±0.3 

— 0.9±0.2 

9.2±0.7 

10.9±1.6 

Vaasa 

0.8±0.2 

-0.6±0.1 

8.8±0.5 

0.4±0.2 

— 0.4±0.1 

8.4±0.5 


Vilhelmina 

— 1.5±0.1 

0.8±0.1 

6.4±0.6 

— 2.1±0.3 

0.0±0.2 

5.4±0.7 

7.3± 1.6 

Virolahti 

0.3±0.3 

— 1.2±0.2 

0.3±0.9 

0.9±0.2 

— 0.8±0.1 

— 1.1±0.4 


Visby 

— 0.8±0.1 

-0.8±0.1 

2.0±0.6 

— 1.3±0.2 

— 1.7±0.1 

— 0.9±0.6 

1.8±1.3 

Vanersborg 

1 

O 

O 

H- 

o 

i — 1 

-0.3±0.1 

2.9± 1.0 

-1.3±0.4 

0.2±0.2 

6.0±1.0 

4.8±2.1 

Ostersund 

— 0.6±0.1 

0.5±0.1 

6.3±0.7 

— 2.3±0.2 

O.liO.l 

7.9±0.6 

7.6±1.7 

Overkalix 

-O.liO.l 

2.1±0.1 

2.9±0.9 

1.8±0.3 

1.7±0.2 

4.5±1.0 

7.7±2.0 


The uncertainties shown are the standard errors (see text). 
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Table 3. Tidal loading effects. 



Effect 

Assumed 

Amplitude,* 

mm 

Example 

Sites 

Calculated 

Admittance 

Total 

Amplitude. 

mm 

(a) 

Snow, soil 

300 

Inland sites 

-0.035 

-10 


moisture 


Coastal sites 

-0.021 

-6 

(b) 

Kattegat! 

±500 

Onsala 

-0.015 

±8 


water level 


Vanersborg* 

-0.005 

±2 

(c) 

Baltic 

±500 

Visby 

-0.020 

±10 


water level 


Oskarshamn 

-0.016 

±8 




Finnish sites 

-0.012 

±6 

(d) 

Great lake 

±10 3 

Karlstad 

-0.004 

±4 


water level 


Jonkoping 

-0.003 

±3 

(e) 

Hydropower lake 

±10 4 

Coastal sites 

-0.001 

±10 


water level 


Other sites 

30%(?) 

±(3-5) 


* Estimated from various sources. 
*Also Hassleholm. 
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Table 4. Pillar measurements. 


Pillar 

Date 


Reference Position* 


North, 

mm 

East, 

mm 

Up, 

mm 

Arjeplog 

93/08/16 

-0.164 

0.208 

458437.694 

Kiruna 

93/08/17 

4574.341 

-4558.030 

2900.038 


94/06/15 

4574.454 

-4558.301 

2900.163 


94/06/15 

4574.249 

-4557.917 

2900.015 


95/06/16 

4574.347 

-4558.288 

2900.452 

Skelleftea N 

93/08/15 

-0.366 

0.434 

58916.931 


95/06/15 

-0.612 

0.249 

58916.984 

Skelleftea S 

93/08/15 

0.817 

1.487 

58941.716 


95/06/15 

1.107 

-0.139 

58942.120 

Hassleholm 

93/06/14 

0.001 

-0.001 

78509.726 


95/05/19 

-0.012 

-0.220 

78509.310 

Jonkoping 

93/06/18 

0.006 

0.039 

227368.024 


94/06/23 

0.028 

-0.070 

227367.941 


94/06/23 

0.002 

-0.151 

227367.891 


95/05/21 

-0.186 

-0.233 

227368.054 

Karlstad 

93/08/12 

-0.200 

0.026 

82786.824 


95/05/23 

0.042 

0.563 

82786.121 

Leksand 

93/08/11 

-0.264 

-0.253 

447573.368 


94/01/18 

-0.278 

0.063 

447572.154 


94/01/18 

-0.276 

0.121 

447572.182 


94/02/06 

-0.327 

0.145 

447572.615 


94/02/06 

-0.335 

0.114 

447572.491 


94/03/08 

-0.428 

0.097 

447572.502 


94/03/08 

-0.431 

0.158 

447572.462 


94/04/15 

-0.576 

0.218 

447572.577 


94/04/15 

-0.590 

0.176 

447572.433 


94/06/14 

-0.421 

-0.336 

447572.718 


94/08/25 

0.252 

-0.357 

447572.786 


94/08/25 

0.263 

-0.304 

447572.767 


95/10/05 

-0.363 

0.146 

447572.943 

Martsbo N 

93/08/05 

0.656 

0.778 

50551.547 


95/10/06 

0.098 

0.570 

50550.558 
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Table 4. (continued) 


Pillar 

Date 


Reference Position* 


North, 

mm 

East, 

mm 

Up, 

mm 

Martsbo S 

93/08/05 

-1708.119 

6.033 

50546.508 


95/10/06 

-1707.842 

5.270 

50545.648 

Lovo 

93/10/28 

-0.537 

-0.499 

56086.280 


95/05/16 

-0.532 

-0.215 

56086.310 

Norrkoping 

93/08/22 

-0.074 

-0.046 

12870.190 


94/06/29 

0.540 

0.481 

12870.061 


94/06/29 

0.425 

0.308 

12869.927 


95/05/17 

0.286 

0.303 

12869.456 

Onsala 

93/08/16 

-216.580 

-171.680 

10039.377 


95/05/20 

-216.127 

-171.758 

10039.511 

Oskarshamn 

93/06/16 

-0.062 

0.017 

119433.397 


95/05/18 

-0.170 

-0.221 

119432.693 

Sundsvall 

93/08/09 

-0.091 

-0.300 

7094.946 


95/06/13 

0.754 

-1.390 

7094.823 

Sveg 

93/08/03 

0.029 

-0.012 

458203.058 


95/06/21 

-0.040 

0.077 

458203.043 

Umea 

93/06/13 

0.000 

0.003 

31531.774 


95/06/14 

-0.544 

-0.665 

31534.045 

Vilhelmina 

93/08/11 

0.043 

-0.063 

420142.000 


95/06/18 

-0.394 

0.951 

420141.126 

Visby 

93/08/12 

-1.866 

1.323 

54864.592 


93/08/12 

-1.514 

1.154 

54863.358 

Vanersborg 

93/09/09 

0.376 

0.605 

134762.710 


95/05/22 

0.092 

0.975 

134762.404 

Ostersund 

93/08/10 

-1.045 

0.285 

458386.234 


95/09/15 

-0.19 

0.009 

458386.585 

Overkalix 

94/06/15 

0.517 

0.413 

200290.683 


94/06/15 

0.474 

0.301 

200290.535 


95/06/16 

0.472 

1.064 

200290.406 

* Heights 

in the Swedish national height system. 

Northing and easting 

in local 


reference system established at each site. 
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Figure 1 . Map showing GPS sites used in the study. Subnetworks are indicated 
by the symbols used: SWEPOS (triangles) and FinnRef (circles). IGS sites are 
indicated by diamonds. The Kiruna IGS site, located close to the Kiruna SWEPOS 
site, is not shown. The Onsala and Metsahovi sites are also IGS sites. 

Figure 2. Time series of site positions for the BIFROST sites and Tromso from the 
Standard Solution. Shown in blue are differences in mm from a series mean for the 
east (left column), north (middle), and radial (right) components of site position. 
The red line shows the model fit to the time series, and the yellow vertical lines 
indicate epochs at which an offset was introduced. 

Figure 3. Comparison between Standard, Edited, and EOF solutions. 

Figure 4. (a) Vertical crustal rates from GIA predictions calculated using a method 
outlined in the text, (b) Simple Gaussian model produced by a least-squares fit to 
the observed GPS vertical crustal rates. 

Figure 5. Comparison between the rates observed from the GPS data and those 
predicted by the Earth/ice model combination discussed in the text for the (a) east, 
(b) north, and (c) radial components of velocity. 

Figure 6. Observed sea-level rates from Baltic tide gauges versus negative crustal 
vertical rates obtained from the Gaussian model of Figure 5b. The mathematical 
relationship between these quantities is given in (2). 

Figure 7. Residuals for radial component for Umea from EOF analysis (a) first 
and (c) second stages, after the recycling of the (b) “common mode” based on the 
residuals for the network of sites. For clarity the error bars are not shown. The 
WRMS residuals are (a) 7 mm; (b) 4 mm, and (c) 3 mm. 
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Figure 8. Eigenvectors in model space. The signal channels through which resid- 
uals at neighbor stations have been recycled are numbered 15-25. The associated 
GPS stations are, in order of their parameter number, Metsahovi, Skelleftea, Uniea, 
Sundsvall, Martsbo, Leksand, Sveg, Ostersund, Vilhelmina, Arjeplog. and Kiruna. 
The residuals contribute to the set of signals from which the (orthonormal) eigen- 
vectors are constructed. Eigenvector number 01 conveys a common mode of the 
recycled noise. Deleting from the model space all eigenvectors that convey the 
residual information except the one that is associated with the largest eigenvalue 
constitutes the extended EOF method used in this paper. 

Figure 9. Correlation between vertical time series from pairs of GPS sites as a 
function of intersite line length. The smaller “branch" at a correlation of ~0.2 
contains correlations with site Metsahovi. This site was used to establish the 
reference frame, so the correlated variations in position might have been removed 
to a greater extent . 

Figure 10. Examples of power spectra for the three components of position for 
six sites. The dotted line indicates the 95% confidence level for white noise er- 
rors whose standard deviations are given by the standard errors for the position 
estimates. 

Figure 11. Time series of the pillar measurements for the Leksand site. The error 
bars are equal to the RMS based on the repeated measurements from other sites 
(see text). The symbols represent the different components of position: square 
with solid line, north; circle with dashed line, east; and triangle with short dashed 
line, up. Each of the time series have been de-meaned, and slightly offset in the 
abscissa for clarity. 



imm. : : 


Kiruna* 


Sodankylla 


. r:r v-jS'ff"- 


Arjeplog* 


Overkalix 


Kuusamo 


65°N 


Vilhelmina* 


Skelleftefi 


Oulo 


Ostersund* 

Sveg A 

Leksand* 


Unjei 


Sun^svaff 


MArtsbo 


Kivetty 
Vaasa • 


Olkiluoto 


Romuvaara 


Joensuu 


Virolahti 


60°N 


i 


Tuorta *Me*aShovi 

^Karlstad ^Lovo 
Vinei^borg ^Norrkoplng 

Bords A ^° nk °P> n ^VIsby 

Onsala A Oskarshamn Riga 


55°N 



10°E 


15°E 20°E 25°E 30°E 


Figure 1 














r 


1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 


1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 




r 


1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 


1994 1995 1996 1997 1998 


1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 

Yilhelmina i r 


■Miu 


1994 1995 1996 1997 1998 


1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 

Virolahti i 



[ 1 1 1 u 1 

1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 

1994 1995 1996 1997 1998 
Figure 2 (continued) 





BIIILkIi 


1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 

Ostersundi i 


PPPi I 


1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 

' 1n — Overkalixi t • 


1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 1994 1995 1996 1997 1998 


Figure 2 (continued) 



Rate from 





66 


A B 



10’E 20°E 30°E 10°E 20'E 30’E 



-3 01 23456789 10 

uplift rate (mm/yr) 


Figure 4 












72 



Figure 10 
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